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ABSTRACT 


The generation and propagation of small amplitude 
acoustic waves in a homogeneous, loss-free, compressible 
im.d 1S studied by the finite element method. A diaphragm 
iemaced in an infinite rigid baffle generates acoustic 
Weves in a semi-infinite fluid region. Steady-state pres- 
Seee O1sStribution is found for a hemispherical region with 
memnaary reilection suppressed through use of a radiation 
femaitton. The computer program developed for the purpose 
Mai ZeS +60—parametric finite elements with curvilinear 
boundaries. incorporated in the program is a versatile 
Mee fsenerator which minimizes the quantity of input data. 
Peeeepvable agreement with analytic results is obtained when 


there are at least four elements per wave-length. 
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1. [RO DUC ae 


ihe fFiela of acoustics fetays an inpOorvane ole im ine 
development of the complex sonar systems used in The United 
scates Navy today. It is the intent of thiS paper to open 
anew avenue to the study of the active sonar system. 
meeemcvion is confined to the transducer and the periodic 
pressure field induced in a fluid medium by harmonic 
@ecililations of the transducer piston or diaphragm. 

mje SUtUCyemamrestricted to the steady-state performance 
Gia transducer mounted in an infinite rigid baifle and 
radiating into a half-space filled with a homogeneous fluid. 
Both the pressure field and the acoustic impedance of the 
transducer are determined. Figure 1 illustrates the 


geometry of the region. 
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Figure 1. Geometry of Region. 











The method used in this study is the two-dimensional 
finite element technique. Since it is manifestly impossible 
wo. represent an infinite regi On wages, Piero e ump) sem 
fanite elements, the stratagem of the anechoic chamber is 
emulated mathematically. This is accompixsmed by using 
radiation (non-reflection) boundary conditions. 

fine diaphragm, being cireular, forms a surtace you 
revolution when deflected. Thus, the instantaneous pres- 
mame) GiStribuvion has rotational symmetry about an axis 
fmeoughn the diaphragm center and normal to the rigid bafile 
(@2)oxis of Figurel). <A finite element assembly representing 
ees cCeion Of imverest may be constructed by dividing the 
eermresponding portion of the r-—-Z plane of Figure 1 into sub- 


feeas aS Shown in Figure 2. 
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Figure 2. Finite Element Grid 











In Figure 2, an individual rectangular element represents 
the cross-section of a toroidal element generated by revolu- 
mon about the 2 axis, The wrovaticwSs ee yimetry permis 
calculation of element properties using two-dimensional 
technique by simply taking the thickness (normal to the 


plane of the figure) proportional to r. 





Ti. INIDE SVEN TOR EA ne hie i 


ee = 


“THE GOVERNING EQUATION 








fee EQUATION OF FLUID MOTION 

The governing partial differential equation for the pres- 
sure p in a homogeneous, inviscid, compressible fluid as 
used in acoustics (e.g., see Ref. [1]”) is 


Vp = 5 p (1) 


© 
where We is Ctnewbeplacian 7Oweraror., «esis the acoustic 
velocity, and superior dots represent differentiation with 


Meee ct tO time. For the purpose of the present study, 


Bemmcdary conditaons take the form 


dp : 
In 7 PYy (2) 


where the coordinate n is measured along the normal to the 


memeaary Surface, p is the fluid density, and Vas iS ee 


memeonent of fluid particle velocity in the n direction. 


B. SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS 

Gne finice element method replaces Eqs. (1) and (2) by 
a set of simultaneous ordinary differential equations with 
independent variable time and dependent variables consisting 
of pressures at a large number of discrete points (nodes) 


@e the region. The region is divided into a number of 


“Numbers in brackets designate references listed on page 68. 








sub-regions called "elements." Within each element it is 
postulated that the pressure p may be expressed in terms of 
the pressure ps at the n nodes of the element. This rela- 


tion is expressed as 
n 
Pp =) Se. (3) 


where the "Shape" (or interpolation) functions N. are 
functions of the spatial coordinates. If the nodes are 
itoecaved at suitably chosen points on the element surface, 
unique sets of shape functions will guarantee the necessary 
Sem@oinuity of pressure p along interelement boundaries as 
Meancveda out by Zienkiewicz [2]. The elements used in this 
study are two-dimensional eight-noded iso-parametric 
elements. Further details are fiven in Appendix C. 

It has been shown by Zienkiewicz and Newton [3] that 
the finite element equation for k nodes, m of which are 


ieaved on the fluid-diaphragm interface, may be written as 
= TPS 
Otel Dito) Lhlipl = = oll 5) (4) 


fere ip} is a k x1 vector of nodal pressures, [Q], [D], and 
[H] are k xk symmetric matrices, {6} is an mx] vector of 
diaphragm displacements, anda high is a kK xm matrix. Elements 
eels), [D], [HH], and Mile are all real constants. Formulas 
jor calculating these matrices are given in Appendix E. fhe 
terms [D]{p} and [L] {8} represent boundary conditions and 


fee! be discussed in following sections. 








C. DIAPHRAGM BOUNDARY CONDITION 

The diephragm displacement 6 ¥s sa prescribed tune memos 
time. It can be represented by choosing a number of points 
(nodes) on the diaphragm and using a set of shape functions | 


N TO interpolate the local displacement as 


NS, (5) 


Or 
I 
les SS 


where the ee are the displacements at the m diaphragm nodes. 

The boundary condition of Eq. (2) is formulated in terms 
erevhne normal component Ve Gteriuid particle velocity. 
the diaphragm surface this must, in the absence of cavita- 
tion, be equal to 6. Using this equality, the fluiad- 


Meo nracm JInvertace condition fives the terin 

~p[L)" {8} 
in the finite element equation (Eq. (4)). The resulting 
vector has nonzero elements only for those pressure nodes 
ineecaved on the diaphragm surface. It is noted that the 


Meuncary condition on the rigid baffle 1s a special case, 


ieecnat the particle velocity Ve 0, which implies eO 


D. RADIATION (NON-REFLECTION) BOUNDARY CONDITION 

A normally incident wave meeting the interface between 
myo media will not be reflected if the specific acoustic 
impedances of the two media are equal. The specific acoustic 


aimpedance Zz is defined to be the ratio of pressure to 


10 





article velocity. For plane) Marnonie ya es Ve aia 


Quantity and the pressure can be given as 
Pea oc (6) 
Bor a wave progressing in the positive n direction. Differ- 


Paviating Eq. (6) with respect to time and using this result 


to eliminate Vv from the boundary condition (Eq. (2)) gives 


gP = _ =p (7) 


on 
Equation (7) then becomes the boundary condition for non- 
reflection of normally incident plane waves. The term 


[D}] {vp} 


masa. (41) results from this criterion. This vector has 
nonzero elements only for those nodes which lie on the non- 


meawecting portion of the boundary. 


JL 








Lil. GROMER Bot she el] 


The initial solution of the problem Using the geometric 
monn t.euration of Figure 2 resulved in Signiticame serror in 
the pressure distribution when the solution was checked 
@eaanst known results. The principal cause of this error 
was believed to be the shape of the non-reflection boundary. 
By the nature of the geometry, the harmonic oscillating 
Giaphragm will create progressive pressure waves which | 
Meoparate outward and strike the boundary of Figure 2 at 
Memrous angles of incidence. The particle velocity, Vi? 
will not always be normal to the non-reflection boundary 
ama the boundary condition required by Eq. (2) will not be 
wae tiecd. Thus, a modification to the representation of 
ges oOoundary to provide substantially normal incidence led 
Menthe feometry illustrated inFigure 3. For this modified 
fmearabion boundary, in three dimensions the surface is 
hemispherical. | 

it order wo represent adequately the circular arc of 
Pie non-reflection boundary, two-dimensional iso-parametric. 
fanite elements are used. The element thickness normal to 
mame lane of the figure is directly proportional to the 
fmedius r. Ilso=parametric elements are further discussed 
iy Appendix C, 

Neve Che wdirstGrtron OF Che Brebasrm treure 3. [he Preaen 


mor Chis distortion and the method of generating the grid 
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Figure 3. Modified Finite Element Region 


Pevtern is discussed in Appendix D. The distorted 
"rectangles" still represent the cross sections of toroidal 


elements. 


iS 








IV. WAVE THEORY AT RADIATION BOUNDARY 


A one-dimensional study of a pulsating sphere acting on 
meerauid medium resulted in a more accurate solution when 
spherical wave theory at the non-reflection boundary was 
mmeroauced. Based upon the success Of these results, a 
further modification was made to the computer program 
(Appendix B) to incorporate spherical wave theory in the 
two-dimensional study. This was accomplished by developing 
time NOon-reflection boundary condition for spherically diverp- 
ing waves as shown below. 

The specific acoustic impedance for spherically diverg- 


ing waves is given in Ref. [1] as 
ae ee (oe) 


Mmimere WwW 1S the circular frequency and r is the radial 
Gistance from the "center." The particle acceleration Vn 


can be written as 


oe = juv,, (9) 


feo allows the boundary condition to be represented by 
oe Ojwv (10) 


or i 


Mrom Eqs. (8) and (10), the boundary condition is 
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Oy } WD 
or C W) (11) 
—7 
ong 
op. =e uf - W 
= P(p + ag (12) 


The imaginary term on the right hand Side or BomeauGl: eu 
mepresented by the term [D]{p} = jw[D]{p} in Ba. (4). Thus, 


the contribution of the real term can be written as 
Cc 
=e) (13) 


Pa@cecan be added to the system of equations given in Eq. (1). 
Hine spnerimecatecOonLribution is only added to the nodes 

fmoeen lie on the non-reflection boundary, thus providing a 

spherical wave approximation to the plane wave theory at 


Mmemnon-retflection boundary. 


ib. 








V. COMPUTER SOLUTION TECHNIQUE 


Assume the diaphragm is oscillating with a harmonic dis- 


placement represented in the complex form by 


fa) = (Nes (11) 


where {6} is a vector of prescribed displacement amplitudes 
at the m nodes of the diaphragm-fluid interface. Elements . 
et i} aAremredime he presetee in che fluid medium vo meecm 


meoresented by 
{p} = {p Jel" (15) 


where UD ae Tomomiex  Veeclor Of pressure alplitudes (conple) 
at the k nodes of the region. 


fiiserting Bas. (14) and (1!) into Eq. (4) gives 


({H] - w°[Q] + jw[D]){p} = w°{B} (16) 


where {B} = [L] {6}. 

In the computer program, the system matrices [H], [9], 
iD), and {B} are assembled from the contributions of each 
Mmm@aavidual finite element in various subroutines. since 
only a small number of nodes are located on the non- 
reflection boundary, the matrix [D] is sparse. In order to 
reduce computer storage requirements, [D], which is symmetric, 
is stored as three vectors: the principal diagonal and the 
two diagonals immediately above. At the same time, the 


mooroximacted spherical contribution ={D] at the 


16 





non-reflection boundary is computed. These quantities are 
Pe a1 andware added dipeeci poke einen eee 
fiitatve Jocations. 

The final subroutine introduces the value of omega and 
multiplies the matrices [Q]), [D], and {B} by the appropriate 
power of omega before the left-hand side is summed to form 
one complex matrix. Thus the system of equations is 


reduced to 
[Si fp tt | aly) 


wamere [S} is a kxk complex matrix and {F} is a column 
Meecvor of order kx lin which only m elements @re nonzero. 

The complex pressure amplitudes in Eq. (17) are found 
meus ing the Gaussian elimination and back substitution 
method. 

The absolute value of the pressure and the phase angle 
mmeemcalculated in the program and printed out in array 
form. Additionally, a system energy balance and the 
Seapharagm acoustic impedance are determined. 

The program provides for eee ens paren SOUL On. DY 
using incremented values of omega so that a range of 


frequencies can be studied. 


1i/ 





V ee ele 


Problems studied were chosen to Verity (ine waccirac mot 
mae method through comparison with available, exact 
Menelytic) solvtions. Results for each problem are 


discussed below. 


meoolem 1. Paston in a Rigid Tube 

me verifiveeme correctness of the computer program, 
finite element solutions were obtained for the problem of 
meeremonicallyeescillating circular piston at one end of 
a fluid-filled, rigid pipe of equal diameter. An infinite 
length was simulated by terminating the finite element 
region with a plane non-reflection boundary at a distance 
mimetic DIStTONFdrameter. This is, of course, a one dimen- 
sional problem with uniform pressure amplitude throughout 
and phase angle is a linear function of distance along 
fe pipe. 

Beveailedmeesvlts are not presented fPoervenmis essentially 
Mmemvial test problem. It was found, however, that there 
was excellent agreement with theory for long wave-lengths. 
Bor wave-lengths less than approximately four times the 
element width (measured parallel to the direction of wave 
propagation), the accuracy deteriorated. This is under- 
standable, since the shape functions employed represent 
the instantaneous pressure profile across an element by a 


parabolic arc. The lower limit on wave-length corresponds, 
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for a given element grid; to anvupper amet, Gh t regency. 
Eaoblem 2. Fisvon in an 21 fangs Geek 

The problem explored most extensively is that of a 
meemonitcally oscillating circular pis von eumrcounced by yam 
miginite rigid baffle and radiating into @ semil—intinitve 
med region. For this problem, the AGrsRO MMF 1 Olde Sowvalelsnene TS 
a hemispherical surface with the center located at the 
PeeswOn center and the radius is four times that of the 
Mmiecon. The element grid.(see Figure 4, Appendix A) contains 
16 elements. It is derived by mapping a 4 x4 square array 
maeom the €,n plane into a quadrant of a circle. 

PWalVGvemipestleS tor thas problem are readily avail— 
pole (e.g., see fl}, [4]). In Table I, the pressure ampli- 
tuces and phase angles along the Z axis, obtained by the 
finite element method (FEM), are compared with exact 
results. Two sets of finite element results are fiven. 
For those designated FEM1, the spherical wave correction at 
@eeeradiation boundary is imeluded. For FEM2, the correc-— 
meen is omitted. Details are given in Table I for two 
circular frequencies, w = .1 and .6 rad/sec. The corres- 
ponding wave-lengths are approximately 63 ft and 10 feet. 

Examination of Table I shows that sal results give 
generally good agreement with exact results. Those for 
FEM2 are considerably poorer, although the Paice 
are smaller for w = .6 rad/sec. This is understandable since 
the W/reconuurtnWerone yon Speci tie sacousyle a2ipeccenccs ls 


relatively less important as w increases (see Fq. (12)). 


Ie 





it is believed that a more refined finite element mesh 
(more elements) should afford significant improvement. 
eance the computer program for this problem required 
Be 0K bytes of core storage, it was Mev Tesciple ce 
mevestigate This without major program revisions. 

im predictins transducer overformance Ghee resne tanya ium 
ieeaction on the piston is needed. This is normally reported 
ieee iorm of piston radiation impedance, which is complex 
eeeeis found by dividing the piston force by the piston 
Memecity. lhe real part of the radiation impedance is 
eeeeled the radiation resistance and the imaginary part is 
@euled the radiation reactance. 

able Ji previcges a comparison with exact values for 
igemradiation resistance and reactance. Results are given 
coe hM] and BEM@ over a range of frequencies. Again it is 


feeeerved that tehe spherical correction significantly improves 


the results. Note that the radiation resistance obtained 
by FEM1 continues to give good results through w = .9 rad/ 
oCce. 


Computer ~olulions of this problem were obtained for a 
range of frequency from w = .1 rad/sec to w = 1.2 rad/sec 
by .1 rad/sec increments. Comparisonsat the upper end of 
fais range are omitted from Tables I and Il because The 
departure from cxact results does not permit useful employ- 
ment of the finite element method without grid refinement. 


meme leve Compucver resulvs for pressure emplitvude, phase 
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dee, lela 
PRESSURE AMPLITUDE WAND. PHASE ANCE SON ears 


me ea re ee a se 
——— —————=— ——————— SSS SO SL ee 


PRESSURE. AMEE Dis PHASF ANGLE 




















2 | EXACT FEM2 PEM PEM2 
Pe | eee 1b/ft* dee deg 
ol 0 02061 | .01792 Ph | eee 
I 2 -00889 | .00676 siete 
al I 00469 | .00376 Mer ees 
dl 6 00329 | .00337 le oo 
i g 00247 | .00331 Peay ee 
£6 0 Pied s S625 ?, 221s 
ts D | epee 273 | LRG. 2 
6 I ese » 186 pe Se S015 
8 6 mel 7 por Oy ee Oa 
6 8 .089 O84 ,2| =260n4 
TABLE IT 


PISTON RADIATION RESISTANCE AND REACTANCE 
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RESISTANCE REACTANCE 












EXACT FEMI FEM? ly, Ge AL WEM 1 
a omecwieD. Secu sec | lbisceumle: cee 
a0 /Sec ft NG Tt ft ft 





A, De eh, Ie 
Le H, Bi Die 
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18 one Oe ie. ‘fe 
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Berit — Finite element results with spherical correction. 
BEM2 —- Finite element results without spherical correction. 


Phase angle given is for pressure relative to piston 
displacement. 


Paid density = 1 eens Piston Radius = 2 ft. 
meoustic Velocity = 1 it/séc. Piston Amplitude = 1 i. 


ea 


a 








angle, and impedance components at several frequencies are 
Biven in Appendix A. 
meoblem 3. Paraboli cmDiaphra gma aie senaseme en lee 

This problem, a modification of the peevtove prop len 
meplaces the piston by a circular diaphragm having a fixed 
edge and having a deflection surface in the form of a 
Paraboloid of revolution. Study of this problem was under- 
taken in an effort to improve agreement in the energy 
balance check (discussed below) INCOPPOraved Gnteww i. 
mieoeran. itt was believed that elimination of fluid 
fewwocity discontinuity at the edge of the piston in problem 
iment resuly in a significant improvement. 

Mnalytic results for this problem are not available in 
Sigoncard acoustic texts, but conventional technicues can 


Bemutilized to obtain the pressure p along the Z axis in 


complex form as 


WS (Oh 











a = oa -j— 
meee 20 | pes |e ofa aS ae [ie a (18) 
P 2 JG oe 
Wa 2 
where 
Sano Sep laccMehumamelivude at diapmmasm cenver 
a = diaphragm diameter 
Pe | i! ~ —_ 2 2 * 
Seen. distances (4 + a) 


For this problem only a single finite element solution 
OeeM3) was obtained. The spherical correction is ineluded. 


Mole [il gives a comparison with exact results (Eq. (18)) 
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TAB LB sf ti 


PRESSURE AMPLITUDE AND PHASE ANGLE ON AXIS 














|| PRESSURE AMPLITUDE | PHASE ANGLE | 
Ww 7 EXACT FEM3 EXACT FEM3 
rad/sec ft Ib/et- lb/ft” deg 
al 0 01427 = 
ml 2 00461 21D 6 
a y 00242 -24.2 
1 6 S ONOMNs C se DS 
al 9 00125 ele i 
“e 0 24934 ~23.8 
6 2 1590 -75.5 
6 M ONE He -~142,3 
é 6 0584 = 00re 
% 9 0436 =278.1 





FEM3 - Finite element results with spherical correction. 
Phase angle is for pressure relative to diaphragm displacement. 
miata density = 1 have DOL 
Meeustic velocity = 1 ft/sec. 
ieeaponragm radius = 2 ft. 

Diaphragm amplitude (center) = 1 ft. 


for pressure amplitudes and phase angles for values of 
H® = .1 and .6 rad/sec. 

Examination of Table JII shows that agreement near the 
mop hragm (small values of Zz) is not as good as FEMI in 
Maple IT. At the radiation boundary however, the FEM3 errors 


are comparable to those for FEM]. 





Bmeiey. Balance Cheek, 


Prior experience, by others, suggested that an energy 
balance check would be a useful fTeauure of @ finite element 
program. Specifically, since solutions represent steady-— 
State behavior, the input power, averaged over a cycle, 
eaeuwld equal the corresponding average value of the output 
Meeer. input occurs at the piston, or diaphragm, and output 
mmmeene radiavion boundary. “Expressions for these quantarvies 
eeees GCerived and incorporated in the computer program. 

A number of program errors in the energy balance calcu- 
ie1ons (Appendix B) prevented valid results until the 
meme! Stage of the investigation. Thus, the values of 
average power in and average power out are not included in 
the computer output (Appendix A). Recent computations, 

Meee propram errors corrected, have shown excellent agree- 
Heme OL these Quantities to at least six significant digits 
imecne lower range of the circular frequency. Even at 
higher frequencies, where solution accuracy (pressure 
amplitudes and phase angles) deteriorates markedly, five 


Significant digit agreement is maintained. 
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VII. CONCLUSIONS 


it has been shown that thes ini Cemetemenc sie vy hecere 1 
provide acceptable results concerning generation and prop- 
agation of acoustic waves. The parabolic iso-parametric 
element provides an economical, accurate means of repre- 
senting regions having curved boundaries. 

ine primgipal imitation ol the method) is evigencs ion 
m@percmpirically deduced results that, for acceptable 
accuracy, a single element must not span more than one- 
quarter wave-length. For wave-lengths of the order of 
Mmeenicvude of the piston radius, the number of nodes 
me ouarea would be around 900. For the present program, 
which already utilizes 320K bytes of core storage for 65 


Mees, this is far beyond capability. 
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APPENDIX A 


COMPUTER “OU Eu 


The rows of pressures and phase angles shown in the 
computer output correspond to the numbering system shown on 
mee tinite element region, Figure 4. 

tne Zeros omimred in Thevarray Of Pressures sate aaa. 
angles do not indicate computed values, but facilitate 


memiving the output in array form. 


Row 
Row 
Row 
Row 
Row 
Row 
Row 


Row 





Row 


Figure 4. Format for-Pressure and Phase Angle Array. 
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The output data are arranged accor ein. Om to eeeree age 


sequence. 


MeMl. Piston in infinite baffle wer sper cece Cle) 
added at the non-reflection boundary. 

Meee, Piston in anfinite baffle without spherwcal correc— 
PiLONweGg@ee ay the non-=rentleccion peolidaw 

Mes. Diaphragm in infinite baffle with spherical correc— 


Tion added at the non-reflection boundary. 


ay 
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APPEND LX C 


TSO-PARAMETRIC ELEMENTS 


The family of iso-parametric elements, recently developed 
by Zienkiewicz, Irons, and others, has brought new versa- 
tility and power to the finite element method. The special 
usefulness of these elements is that they may be shaped to 
Pemiorm tO curvilinear boundaries. Attention is confined 
here to two-dimensional parabolic elements. A full descrip-- 
muon Of the entire family is given in Ref. [2]. 

Cons (ctw cecronpuiar Clement ot Figumes>. The 
element has four corner nodes and four midside nodes, 
mee red as Shown. Points in the interior, or on the 
boundary, may be located in the local (&,n) coordinate 
system. The dependent variable, pressure p in the present 
application, is defined at any point by Eq. (3), repeated 


here 


N.p. 6 3u) 
where 


N, = pl+ nn, ) (1+ 8, ) (nn, + 66, -1) 


mer corner node i, and 


oe epee 


eal 





€=1 





eure 5. Hilement Configuration in Local €.n Coordinates 


for mid-side node i, 


CE, ons) are coordinates of node i, 


and Py is the pressure at node i. 
Note that each shape function N. Nase tie Va puemuniI ty. au siode 
eed AS Zero at all other nodes. 
A curvilinear element, such as ener Figure 6, may 
be produced by mapping the rectangle of Figure 5 into the 


Mm) Diane by the relations 
N.Y; Gt) 


where (x, sy) Eo runec eorleoiameceet dinates eOmenoge 1. 
The shape of the element in the x,y plane is uniquely 


determined by specifying the coordinates Cm for the 
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di 


Figure 6. Curvilinear Elements Mapped 
im cCarvestan Coordimares. 


~ *. 


Bent nodes. The edges of this curvilinear element are 
Merabvolic. 

Adjacent elements "share" the same three nodes at their 
snterface. When mapped to cartesian coordinates, these two 
adjacent elements then have the sane curvilinear shape and 
the same x and y coordinates at nodes. Tidewowov1ees Come 
tinuity in the pressure distribution across element 
boundaries. | 

The region developed for the present stucy necessitates 
the use of iso-parametric elements to provide mie Cl ret lens 
non-reflection boundary. This boundary is constructed 
using parabolic arc segments which are mapped by use of the 
shape functions to "fit" a quarter circle powdery ln Cale 


meSian coordinates. 
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APPEND» 


FINITE ELEMENT MESH GENERATOR 


Two subroutines, included in the compurversorocram of 
Moygendix B, should prove useful in @ Variety or finite 
eaeement problems. oe eo ohee these subroutines subdivide 
the region, systematically assign element numbers and 
meee numbers, calculate nodal coordinates; and produce 
a computer plot of the region to be studied. AS an example. 
consider the region of Figure 7. 

The region is taken to be one large iso-parametric 
element (a mapping of thie rectangle of Figure 5). Its shape 


is prescribed by specifying the (x,y) coordinates of the 





Figure 7. Basic Geometry of Finite Element Region. 
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eight nodes. Subdivision of the element can then be 
accomplished at predetermined Vadwecme! so eancemeo, emia 
Equation (19). 

Subroutine NODE (Appendix B) requires two input 
parameters, the number of divisions desired in the x direc-— 
Fon and the number of divisions desired in the y direction. 
These parameters are Psd to calculate the increment values 
of € and n for subdividing the element into smaller 
elements. The subroutine starts numbering global nodes at 
the upper left hand corner (&=-1, n=1) then adds Cena 
values of € and assigns global node numbers from left to 
iment on the boundary n=l. Similarly, aitter incrementing 
we the second row of giobal node numbers is assigned, etc, 
Global node numbers are assigned to all corner nodes, mid- 
side nodes and the center of each subdivided element (9 
nodes per element). The center node provides unique 
imenbification for the element and simplifies the printing 
of computer output in a readable form. Subsequently all 
nodes are renumbered, omitting the center nodes. 

The cartesian coordinates corresponding to the global 
node numbers are calculated by a similar process in the sub- 
routine XYFORM and are used in the program to plot the 
finite element grid. Figure 8 illustrates a computer plot 
of the finite element mesh generated for a quarter circular 
region (Figure 7) with an eight foot radius. The number 
eo Ssubdivietonssi1m DOvn thew and vy drreeiron 1s equals le 


four. The plotting routine produces straight segments 










1. Ss ee es ee os Oe 
e, ees 


eee Cae ae 


Figure 8. Finite Element Mesh Generated 
by Subroutine XYFORM. 


between nodes. In this instance, the element boundaries 
eee Gill parabolic arcs except for those lying on the x and 


y axes, 
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APPENDIX E 


MATRIX ELEMENT FORMULAS 


Individual elements of the matrices [Q], [DJ], [L], and 
[H] of the system of ordinary differential equations (Eq. (4)) 


ouee Siven by Ref. [3] as 





1 

q.. = —~ f N.N.dAR 

ae) of Re t iy 

ae | Pf ieecls 

at Ca 1 J 
£ . =f N_N.AS 

65) et J 

aN, 9N, aN, 9N, N, ON, 

ng 4 lax ae tay ay te ae | 


Where Re and Se are the element region and external boundary 
respectively. 

tne shape functions N.,; N, are specified in terms of 
local coordinates and it is advantageous COmperrorn the 
integrations in these coordinates. Gaussian quadrature is 


employed. Details of the procedure are given in Ref. [2]. 
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